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Abstract 



We use a modified Shan-Chen, noiseless lattice-BGK model for binary immiscible, incom- 
pressible, athermal fluids in three dimensions to simulate the coarsening of domains following 
a deep quench below the spinodal point from a symmetric and homogeneous mixture into 
a two-phase configuration. The model is derivable from a continuous-time Boltzmann-BGK 
equation in the presence of an intercomponent body force. We find the average domain size 
growing with time as f , where 7 increases in the range 0.545 ± 0.014 < 7 < 0.717 ± 0.002, 
consistent with a crossover between diffusive t 1 / 3, and hydrodynamic viscous, t 1 ' , behaviour. 
We find good collapse onto a single scaling function, yet the domain growth exponents differ 
from others' works' for similar values of the unique characteristic length L and time T 
that can be constructed out of the fluid's parameters. This rebuts claims of universality for 
the dynamical scaling hypothesis. For Re = 2.7 and small wavenumbers, q, we also find a 
q 2 <-> q 4 crossover in the scaled structure function, which disappears when the dynamical 
scaling reasonably improves at later stages (Re = 37). This excludes noise as the cause for 
a q 2 behaviour, as analytically derived from Yeung and proposed by Appert et al. and Love 
et al. on the basis of their lattice-gas simulations. We also observe exponential temporal 
growth of the structure function during the initial stages of the dynamics and for wavenum- 
bers less than a threshold value, in accordance with the diffusive Cahn-Hilliard Model B. 
However, this exponential growth is also present in regimes proscribed by that model. There 
is no evidence that regions of parameter space for which the scheme is numerically stable 
become unstable as the simulations proceed, in agreement with finite-difference relaxational 
models and in contradistinction with an unconditionally unstable lattice-BGK free-energy 
model previously reported. Those numerical instabilities that do arise in this model are the 
result of large intercomponent forces which turn the equilibrium distribution negative. 
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1 Introduction 



Homogeneous binary fluid mixtures phase segregate into two phases with different compositions 
when quenched into thermodynamically unstable regions of their phase diagram, a process also 
called spinodal decomposition. This is achieved by lowering the temperature well below the so 
called spinodal temperature. For incompressible, 50:50 mixtures, also called critical or symmetric 
mixtures, these phases form interconnected domains which at late times produce a bicontinuous 
structure with sharp, well developed interfaces. For asymmetric mixtures (phases with different 
densities) there is a phase transition at early times from an interpenetrating structure termed as 
"bicontinuous" to the so-called "droplet phase" , which in turn undergoes subsequent coarsening 
via coalescence pQ. The composition of a binary immiscible fluid is one of the variables affect- 
ing its dynamics. Fields where spinodal decomposition is of industrial relevance comprise the 
metallurgical, oil, food, paints and coatings industries. Polymer blends and gels immersed in a 
solvent are also potentially important applications where phase separation occurs and needs to 
be controlled [13 CHI 

Spinodal decomposition has been extensively studied by experimental [2], analytical Oil], 
and numerical El EZl IHl El EHl [Til [T2l EEl EH] approaches. The fact that it entails a variety 
of mechanisms that can act concurrently and at different length and time scales has made it a 
testbed for complex fluid simulation methods. Among the latter are hydrodynamic lattice gases 
[46] . the lattice Boltzmann equation and dissipative particle dynamics [47] . 

Despite all the interest attracted by the subject, how the mechanisms responsible for do- 
main separation act remains on unsettled grounds. In particular, the dynamics of the late 
time, true asymptotic growth is unclear. Also, the dynamical scale invariance hypothesis (to 
be explained later on in this paper), which is treated almost as canonical by analytical and 
numerical approaches to solving the continuum, local-thermodynamic Cahn-Hilliard equations, 
has experimentally been proven to fail at least under certain conditions [29] . 

Numerical studies on spinodal decomposition include methods at the macroscopic scale, 
based on the numerical solution of either the Navier-Stokes [501 153] or the Cahn-Hilliard equa- 
tions the mesoscopic scale, where lattice-Boltzmann (LB) methods [ElE^i lattice 
gases 013, dissipative particle dynamics [IS], and Ising [35] approaches are examples, and the 
microscopic scale, with classical molecular dynamics [12] . 

Fluid dynamical methods in the mesoscopic scale came to light as a way to grasp the relevant 
thermohydrodynamical behaviour with as little computational effort as possible. This is achieved 
by evolving a microworld in which the usual vast number of molecular degrees of freedom and 
characterisation have been drastically reduced, based on the fact that, away enough from critical 
points, a fluid's macrostate is pretty much insensitive to many of its microscopic properties. Some 
regard the Cahn-Hilliard equations to be within the mesoscopic scale. They derive from the 
van der Waals' formulation of quasilocal thermodynamics [S3) extended by Cahn and Hilliard 
[51] . and aim at solving a Langevin-like diffusion equation for the conserved order parameter. 
This equation involves a chemical potential derived from a phenomenological, Ginzburg-Landau 
expansion for the free energy, and leads to phase segregation if the temperature is below a 
critical value. The scheme commonly used for the study of phase segregation in immiscible 
fluids is termed Cahn-Hilliard Model H [33] : hydrodynamics is included by introducing mass 
currents, which couple the diffusion equation with the Navier-Stokes equation. Thermal effects 
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are sometimes included in the dynamics by the addition of a noise term satisfying a fluctuation- 
dissipation theorem. 

Cahn-Hilliard equations have been applied to model the segregation dynamics of deep and 
sudden thermal quenches of fluid mixtures. Such quenches are usually chosen to be sudden to 
avoid thermal noise effects and set up an initial condition that quickly leads to a state of steep 
domain walls and where diffusion is negligible compared with hydrodynamic effects, thus leaving 
the conditions that the dynamical scaling hypothesis requires. However, local-equilibrium cannot 
be guaranteed for a mixture undergoing a sudden quench, which puts the existence of a free 
energy and the equilibrium states modelled by it on rather shaky grounds. 

The lattice-Boltzmann method we use in this work is the Shan-Chen lattice-BGK scheme 
for binary immiscible and incompressible fluid flow The equilibrium state for each pure 

fluid is chosen to be a local isothermal Maxwellian, and Shan and Chen's contribution to the 
lattice-BGK scheme comes through the phase separation prescription. This is incorporated via 
intercomponent repulsive mean-field forces between elements of fluids (meant to be at a meso- 
scopic scale) which alter the local equilibrium, and not through a local equilibrium reproducing a 
chemical potential derived from a free-energy functional. The Shan-Chen method has been used 
by Martys and Douglas [31] to qualitatively simulate spinodal decomposition for critical and 
off-critical quenches in 3D. There have been recent quantitative studies in 2D using this method 
for critical spinodal decomposition |22l I3()| . An early study on critical 2D and 3D spinodal 
decomposition was put forward by Alexander et al. j35j using the lattice-Boltzmann method 
proposed by Gunstensen et al |3l)] . 

Lattice-BGK methods based on a Ginzburg-Landau free-energy functional |15j achieve mul- 
tiphase behaviour by using two separate distribution functions: one for the mass density and 
one for the order parameter, this being defined as the difference between the phases' densities. 
Higher-order velocity moments of these distributions are imposed to coincide with thermome- 
chanical quantities obtained from the free energy. The term "top-down" is used in the literature 
to address this type of approach, whereas we shall use "bottom-up" in the remainder to signify 
fully mesoscopic methods. Some criticisms of top-down approaches 59 include their frequent 
lack of Galilean invariance (although Inamuro presented a model that does exhibit this 
property), and their phenomenological character. Studies of spinodal decomposition using these 
methods are described in the works of Wagner et al. ^3 , Kendon et al. 9 , and Cates et 
al. [2S1- 

Numerical instabilities are a great cause for concern in lattice-Boltzmann methods, a study 
of which will be addressed for the lattice-BGK method we employ in this work. Their sources are 
two-fold: (a) the finite-difference, discrete- velocity scheme used to solve the BGK-Boltzmann 
equation prevents the existence of an H-theorem, and (b) the approximations used for the 
equilibrium distribution do not guarantee its positivity, and hence that of the nonequilibrium 
distribution. Linear stability analyses have been applied to the lattice-BGK model by Ster- 
ling [35], and in more detail by Lallemand and Luo 00] comparing a lattice-BGK model to a 
generalised LB model with a different relaxation time for each physical flux. Qian et al. |41| 
gave conditions for the Mach number and the shear viscosity such that the lattice-BGK scheme 
produces positive mass densities. New approaches to unconditionally stable lattice-Boltzmann 
models have recently appeared too |42[ I43| 1441 I45j . They prove the existence of functionals 
satisfying an H-theorem. 
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Our objective in this work is to present a bottom-up lattice-BGK method for the study of 
scaling laws in the spinodal decomposition of critical fluid mixtures in three dimensions. This 
method has certain advantages over lattice-BGK methods based on a free-energy functional, 
namely, a smaller number of free parameters to tune, Galilean invariance guaranteed, and a 
simpler equilibrium distribution. Moreover, it refuses to inject macroscopic information into 
the mesoscopic dynamics as the top-down methods do, on the grounds that for lattice-BGK 
methods there is no H-theorem available that guarantees an unconditional approach to a given 
equilibrium. Indeed, in the context of general complex fluid applications, an expression for the 
free energy itself may be unknown, and/or its validity be questioned for regimes far enough from 
local equilibrium, making a top down approach not even viable. 

The remainder of the paper is structured as follows. In Section [2 we discuss the dynamical 
scaling hypothesis, which asserts that after the quench all length scales in the mixture share 
the same growth law with time. The modified Shan-Chen model we use is explained in Section 
El Section 0] introduces the method we use to measure surface tension, while in Section El we 
describe the simulations performed and the growth laws and scaling functions drawn from them 
which allow to test the validity of the dynamical scaling hypothesis. Finally, we present our 
Conclusions in Section [7| 

2 Spinodal decomposition 

After domain walls have achieved their thinnest configuration via diffusion, the time evolution 
of the bicontinuous structure that is produced in the phase segregation process that symmetric 
mixtures undergo presents geometrical self-similarity to the initial stages of such a process when 
the structure is zoomed in at increasing magnification. This leads us to the dynamical scaling 
hypothesis, which states that at late times, when diffusive effects have died out, there is a 
unique characteristic length scale L which grows with time such that the geometrical structure 
of domains is (in a statistical sense) independent of time when lengths are scaled by L |19j . This 
amounts to saying that all length scales have the same time evolution. Such a characteristic 
length scale must be universal for all fluids with the same shear viscosity, rj, density, p, and 
surface tension a, provided that no mechanisms are involved in their late stage growth other 
than viscous dissipation, fluid inertia and capillary forces, respectively. This is so because, as 
we shall see later on, only one length scale can be constructed out of the fluid's parameters rj, p, 
and a, these being the only ones present in a hydrodynamic description of the mixture via the 
Navier-Stokes equations. 

The characteristic length scale is usually measured by looking at the first zero crossing of 
the equal-time pair-correlation function of the order parameter fluctuations ^Z] , 



where, on the lattice, () = X)x?/^> ^ i s the spatial volume, <; is the volume of the lattice's 
unit cell (hence V/s is the number of nodes in the lattice), T is the time parameter in time 
steps, r and x are spatial vectors, and (f)' = <f> — (4>) are the order parameter fluctuations, where 
0(x) = p R (x) — /0 B (x) is the order parameter for our binary fluid (say, a mixture of red (R) 
and blue (B) phases). The units of C(r,T) are squared mass density. In the remainder, "lattice 
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units" will mean unity for the mass, length, and time units, respectively, in an arbitrary unit 
system. The Fourier transform of C(r,T), called the structure function, is 



5(k,r) = ^(r) 
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The units for the structure function are the same as those for the correlation function, and 4>'^ 
is the Fourier transform of the fluctuations. Function ((2J) is volume-normalised, and gives no 
power spectrum for infinite lengths, i.e. 

^J2S(k,T) = l, S(k = 0,T) = 0. 
k 

where is the unit cell volume in reciprocal space, and V = (2tt / 'L) 3 V '/ 'q = (2-7r) 3 /^ is the 
reciprocal space volume; in fact, s'/V = s/V. Although Eqs. Q and (J2J) are numerically 
equivalent, the intensity of X-ray and neutron scattering is directly proportional to the structure 
function, which is hence easily measurable; it is thus this quantity that we prefer to use to 
measure the system's characteristic length scales. 

We define the (time-dependent) characteristic size L of the domains as 

in lattice units, where k\(T) is the first moment (mean), 

fcl ^ E vS ' (4) 
of the spherically averaged structure function S(k,T), defined by 

5(fc,r)^ E ^ (k ; T) , (5) 

2^k 1 

where k indicates the set of wave vectors contained in a spherical shell of thickness one (in 
reciprocal-space lattice units) centered around k, i.e. such that n — \< -k < n + |, n being 
an integer, k is the modulus of k which is smaller than the Nyquist critical frequency k c = n to 
prevent aliasing. In the limit of short distances and large momenta, scaling arguments lead ^5] 
to the relation 

valid for kL 3> 1, also known as Porod's law, where D is the spatial dimension. Short distances 
here means £ <C r <C L, where £ is the interface thickness. 

Other measures have also been used for the system's characteristic length scale, namely, the 
position of the structure function's maximum, and the structure function's second moment, &2 
|17j . We chose to use the first moment k\ as it is the simplest quantity among the aforementioned. 
Appert et al. [HI] found that the structure function's maximum's wavenumber provided a length 
evolving similarly, although in a noisier fashion, to that derived from the first moment. 
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Mathematically, the dynamical scaling hypothesis can be written as 

C(r,T)=f(r/L), (7) 

or 

S(k,T) = L D g(kL), (8) 

where L = L(T) is a function of time, and g is the Fourier transform of /, both of which are the 
same for any late time slice. 

Using methods introduced by Kendon et al. 0, there are unique length and time units 
that can be defined from the fluid's density, shear viscosity and surface tension, p, r/, and a, 
respectively, as follows: 



£ = — , T = -^. (9) 

pa pa A 

We can think of these as a wavelength and a period associated with the system's fluctuations, 
respectively, although they do not necessarily have to refer to actual fluctuation averages. We 
can define the dimensionless variables: 



I = L/L , t = (T - T int )/T , (10) 

which serve to express the universal character of the dynamical scaling hypothesis. Parameter 
Ti n t is an offset that allows one to account for early time diffusional transients and lattice effects. 
Due to the finite resolution of the lattice the initial condition is not an infinitely fine-grained 
thorough mixture ((f) = 0) but there is a non-negligible domain size measured at time T = 0. 
We have then to specify a time origin prior to T = 0, corresponding to a fictitious zero domain 
size. 

For a critical binary immiscible mixture in three dimensions, scaling arguments applied to 
the terms of the Cahn-Hilliard Model-H equations show that Eq. (JJJ holds in the asymptotic 
limit ^H]) or ) equivalently, that 



/oct 7 , (11) 

where 7 = 1 and 7 = 2/3 for the cases when hydrodynamic viscosity and inertia dominate the 
dynamics, respectively. From the Cahn-Hilliard Model B, which is a Langevin diffusion equation 
without noise conserving the order parameter an exponent of 7 = 1/3 is derived, identical to 
that obtained from the Lifshitz-Slyozov theory for the growth of a minority phase whose volume 
fraction is negligible, and is expected to appear at diffusive stages, before hydrodynamics kicks 
in. Scaling theories do not give any prediction for the crossovers' positions in time other than 
that they are "of order unity" |26j . 

Using a free-energy based, lattice-BGK method, Cates et al. [22] reached the viscous regime 
(I oc t) for Lq w 5.9 and Re < 0.1, and the inertial regime (/ oc t 2 / 3 ) for Lq « 0.0003 and 
Re < 350. The Reynolds number is defined in this domain-coarsening context as 
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(12) 



where v is the kinematic viscosity of the fluid mixture, as defined in the next section, and I is 
the time derivative dl/dt. 

There is also experimental 29\ and 2D simulation JOJ evidence of breakdown of scale invari- 
ance in symmetric binary immiscible quenches. In those experiments, the breakdown of scale 
invariance occurs [22] for symmetric binary mixtures in confined geometries under the influence 
of wetting, and a universality has been reported to hold. The process consists of a hydrodynamic 
coarsening occurring faster than mass diffusion, leaving the system with macroscopic domains 
whose concentrations are near to but not at the coexisting equilibrium ones. Metastability or 
instability of the domains then causes a secondary phase separation to kick in via diffusion. 
Scale invariance and self-similarity have also been recently found to break down for viscoelastic 
binary fluid mixtures jHJ . Finally, there is simulation evidence of breakdown of scale invariance 
coming from free-energy based, lattice-BGK simulations in 2D. The rationale for this is the co- 
existence of competing mechanisms at all times in the mixture: diffusion, hydrodynamic modes, 
and surface tension, giving rise to length scales with different growth exponents [TO] . 

3 Our lattice-Boltzmann model 

Initially introduced as a coarse grained version of the lattice-gas automaton method for fluid 
flow simulation, the lattice-Boltzmann model can also be interpreted as a finite difference solver 
for the Bhatnagar-Gross-Krook (BGK) approximation to the Boltzmann transport equation 
[THj . From lattice gases it inherits a particulate, mesoscopic character, as their particles can be 
assimilated to any physical size which is negligible at a hydrodynamic scale; moreover, unlike 
lattice-gas automata, no fluctuations are present within the scheme |32| . From the simplicity 
of the Boltzmann-BGK collision term the LB method gains algorithmic efficiency in simulating 
fluid flow over solving the incompressible Navier-Stokes equations. When extended to multiphase 
flows, these features are especially valuable in looking at the complicated domain interfaces that 
arise in the coarsening of binary mixtures. 

The method we use is a modification of the multicomponent, immiscible fluid LB scheme of 
Shan and Chen [Tl], which will be explained in detail in Subsection 13.21 The Shan-Chen LB 
model employs an expansion in Mach number of a Maxwellian equilibrium distribution. Phase- 
segregating interactions are introduced by means of a self-consistently generated mean-field 
force between particles. The inclusion of this force gives rise to a non-ideal gas equation of state 
through the Navier-Stokes equation, which is reproduced via the usual multiscale Chapman- 
Enskog [20] or moment (Grad) |21j expansion of the distribution function. No thermohydrody- 
namic behaviour is imposed on the equilibrium distribution, as aforementioned free-energy based, 
lattice-BGK methods do l£j , partly because none of the lattice-Boltzmann implementations re- 
ported in the literature so far exhibit an H-theorem ensuring the existence of an asymptote 
towards a prescribed equilibrium, and partly because a purely mesoscopic, mean-field approach 
is preferred here. The coefficients of the equilibrium distribution expansion are determined by 
the conservation of mass and momentum, the property that Galilean invariance holds, and an 
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isotropic pressure tensor. 

In this work we employ a pseudo four-dimensional lattice, which is the projection onto 3D 
of the D4Q24 face-centered hypercubic (FCHC), single-speed lattice, where the notation implies 
the spatial dimension (4) and the number of vectors linking a site to its nearest neighbours 
(24). The FCHC lattice guarantees isotropic behaviour for the macrosopic momentum balance 
equation |2*lj . 

In the following subsections we introduce our modified Shan-Chen model, first by looking 
at a non- interacting mixture of gases, and second including the mean- field force term which 
gives rise to a non-ideal gas equation of state. Then, we modify the collision term such that 
the Shan-Chen scheme is consistent with that derived from a Boltzmann-BGK equation in the 
presence of a force. 

3.1 Mixture of ideal gases 

The finite-difference, finite-velocity fully-Lagrangian jSH] scheme for the numerical solution of 
the multicomponent Boltzmann equation, 

n£(x + c fc ,i + l)-ng(x,t) = n£, (13) 

governs the time evolution of the A;th velocity's particle number density for the fluid species 
a in a non-interacting mixture of gases. The lattice-BGK collision term is 

aiM^- tM-f*™ , ( u) 

where the time increment and lattice spacing are both unity, is one of the 24 discrete velocity 
vectors plus one null velocity, x is a point of the underlying Bravais lattice, and a = R, B (e.g. 
oil (R) or water (B)). The parameter r Q defines a single relaxation rate towards equilibrium for 
component a. The function n^ eq ^(x, t) is the discretisation of a third-order expansion in Mach 
number of a local Maxwellian [22] , 



a(eq) 



x. /) - -o-'a .//'(x. /) l + ^c fc -u+^|(c fc -u) 2 -^|M 2 + ^|(c fe -u) 3 -^|U 2 (c fc -u) . (15) 



where u>fc are the coefficients resulting from the velocity space discretisation, and c s is the speed 
of sound, both of which are determined by the choice of the lattice. For the projected-D4Q24 
lattice we use, the speed of sound is c s = and ujj- = 1/3 for the speed Ck = and 1/36 

for speeds = l,y/2 [IJ. (The projection from 4D to 3D puts an additional speed into play, 
V2.) InEq. JTS}, u is the macroscopic velocity of the mixture, through which the collision term 
couples the different velocities c^, and is a function of x and t. Also, n Q (x, t) is the local particle 
density for the a-th component, defined as J2k n t( x ^)- 

The judicious choice of the coefficients in the expansion of the equilibrium distribution l)15|) 
allows for mass and momentum to be conserved, 

$>£ = 0, £>«E c * n fc=°- (16) 
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Momentum conservation requires the fluid's macroscopic velocity u to be denned in terms of 
the macroscopic velocity u a for component a, 

n a (x,t)u^]Tn£(x,t)c fc , (17) 
k 

as the solution of the three equations: 

H(u) = v, (18) 

where 

Sj(u) = (2 — 3n 2 )uj + 3uf + 3iijU 2 +1 + 3ujU 2 +2 , (19) 
with the cartesian index i ranging in the imod3 set, and v being defined as the special average: 

"E^/E^- <») 

a a 

Based on previous experience with lower orders, our choice of a third-order Taylor expansion 
in Mach number for the Maxwellian equilibrium distribution is an attempt to improve the 
approximation for velocities which, within the incompressibility limit, are large enough to make 
either the distribution funtion become negative or the error in the expansion too large. 



3.2 Mixture of interacting, non-ideal gases 

In order to deal with non-ideal gases, in particular, fluid mixtures whose volume elements interact 
among themselves, each fluid is forced to relax to a local equilibrium which is modified by the 
presence of its surrounding volume elements. The mean- field force density felt by phase a at 
site x and time t from its surroundings is defined as: 

F «( Xjt ) = -^(x,t)^ fe ^r(x',t)(x'-x) (21) 

a x ' 

where g a a (> for immiscible fluids) is a coupling matrix whose non-diagonal elements control 
interfacial tension, and ip a is the so-called effective mass, which serves as a functional param- 
eter and can have a general form for modelling various types of fluids. For simplicity in our 
implementation, we have chosen ip a (x,t) = n Q (x, t) [22] an d on by allowed nearest-neighbour 
interactions, x' = x + c^. Other choices for if) have also been made |14| . 

Shan and Chen |14j incorporated the above force term in the collision substep of the LB 
dynamics by adding the increment 

Au a = — r a (22) 

to the velocity u which enters the second-order expansion of the equilibrium distribution func- 
tion. We perform the same procedure for our third-order expansion (|15|) . obtaining additional 
terms: 
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n 



a(eq) 



u + Au c 



a(eq) / \ 
«fc ( u ) 



Cfc — u (2cfc ■ u — u 2 



+ 



2c! 



"Cfc 



a" -a" (c fc -a") 2 (c fc • u)(c fc • a") 2 

o a ■ a 



+ —u; k n a (c k -8L a T a ) 3 



,a\2 



(23) 



where a a = F a /p a . 

Luo [HHl and Martys et al. [6j| expanded both the velocity space gradient in the BGK- 
Boltzmann equation force term: 



a • Vt n 



(24) 



and the equilibrium distribution in Hermite polynomials in the lattice velocities. Then they 
rearranged the acceleration a such that it explicitly modifies the macroscopic velocity in the 
equilibrium distribution, leaving a term linear in a. If only linear terms were to appear in Eq. 
(j23|) . the Shan-Chen prescription for an interparticle force would then coincide with the way it 
is included in the continuum BGK-Boltzmann equation, as pointed out by Luo and Martys et 
al.. To this end, following Nekovee et al. (22], we simply drop from Eq. (j23j) any term nonlinear 
in the acceleration a. We thus obtain a modified Shan-Chen collision term, which is why our 
model is termed modified Shan-Chen. The modified Shan-Chen collision term is: 



where 



and 



a I 



o \fiaa^k Saa^l) T L,ace a ^k 



(25) 



(26) 



Co 



(27) 



where we have made use of the condition ()29|) below. The second term arising in (|25|) ac- 
counts for interparticle interactions other than the binary collisions implicit in the Boltzmann 
collision term, £1 This includes a collision operator resulting from mean-field interac- 
tions among different fluid components which gives rise to phase separation for immiscible 
multicomponent systems. 

The inclusion of a mean-field force in the Shan-Chen model leads to the breakdown of the 
local momentum conservation that holds for noninteracting ideal gases, cf. Subsection (|3.1|) . 
However, the forces felt by neighbouring portions of fluid follow an action-reaction mechanism 
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that leads to global momentum conservation (i.e. over the whole lattice). This was numerically 
confirmed for our third-order-equilibrium, modified Shan-Chen model too. 

It can be shown that the condition for momentum conservation in the absence of interac- 
tions, Eq. qi8|). leads to that needed when using a second-order expansion of the equilibrium 
distribution, namely 

u = v , (28) 
only in the limit of creeping flows to second order, i.e. 

u 2 ^0. (29) 

We therefore implemented the computation of the velocity according to (|28|) rather than (|18f) . 
The condition (|29j) is satisfied, as global momentum would not be conserved otherwise. In 
addition, we confirmed in our simulations that the fluid velocity was kept under 28% of the 
speed of sound by 67% of the lattice nodes. This means squared Mach numbers under 0.08. 
This purports to show that the expansion to third order, implemented in this model to extend 
the parameter space for which the equilibrium distribution remains positive, for momentum 
conservation at least adds very little. 

In our LB model, the kinematic viscosity of the mixture is given by 

f = - = c s 2 ( XaTa - o) ( 30 ) 

P a. ~ 

where c~ 2 = 3 for our lattice, T a is the relaxation time of the ath component and x a is its mass 
concentration defined as p a /p EU- For a re gi° n of pure ath component, 

, = i(r-i) (3D 
which also holds for our case of a 50:50 mixture, for which 

r^2x a = T (32) 

a 



2>aT a = 

a 

since all relaxation times are the same. 



4 The surface tension 

The surface tension a arises as an emergent effect due to inter component interactions. It is 
calculated by measuring the components of the pressure tensor P = {Pij} across a planar 
interface perpendicular to the z-axis through the formula 



a 



Pxx [Zj 



dz 



(33) 



where P« is the flux of the ith component of the momentum across a surface perpendicular to 
the jth. cartesian axis. This pressure tensor, consistent with the force 1)21(1. is 
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P W = EEp£( x ) c ^ (34) 

a k 

+ 7 E^E [fWf^l + <Hx)V Q (x')] (x - x')(x - x') , 



a, a 



with x' = x + cj- in this study. This leads to the same expression for the scalar pressure as 
that in the momentum balance equation obtained by multiplying the LB equation ()13|) using 
the collision term (|25|) by and summing over k. Here, /0^(x) is the mass density of species 
a with velocity at the site x. Eq. (34) contains a kinetic term due to the free streaming of 
particles corresponding to an ideal gas contribution, plus a potential or virial term due to the 
momentum transfer among particles of equal and distinct colour, through the interparticle force 

As previously noted, the surface tension in the modified Shan-Chen model is an emergent, 
hence not directly tunable quantity, in contradistinction to the situation with free-energy based 
lattice-Boltzmann models. It depends on the density p, the coupling g and the relaxation time 
r, and has to be determined by simulation. We computed its dependence on these parameters 
to be as follows 

da da da , . 

T P >0 - T 9 >0 - Tt <0 ' (35) 

This behaviour is useful when steering through the parameter space in search of specific values 
of Lq and Tq. Numerical results on the surface tension are reported in the next Section. 



5 Simulations 

We restrict ourselves to critical (50:50) mixtures, which are the type of configurations leading to 
a spinodal decomposition process as opposed to nucleation. Experimentally, spinodal decompo- 
sition is characterised by long-wavelength, infinitesimal density perturbations which are unstable 
after the quench, hence favouring the segregation, whereas nucleation generally presents short 
wavelength, finite perturbations, and metastability is not uncommon. Nucleation is hence a 
more complex phenomenon which is usually considered after an initial study in spinodal decom- 
position has been performed. 

We aim at reproducing the early time diffusive and later time viscous and inertial regimes 
predicted by carrying out scaling analyses on the Cahn-Hilliard Model-H equations |17l I19j . 
Growth laws predicted for those are I oc i 1 / 3 , I oc t and I oc i 2//3 , respectively. Under the 
assumptions of the dynamical scaling hypothesis made in the introduction, those regimes are 
uniquely characterised by the length and time 

L = —^—(r- l -)\ T = —^ r2(r-l)\ (36) 

obtained by inserting (|3*T|) into ©. 
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Having in mind keeping simulation time at a minimum, the values of p, r and g must be 
such as to allow the fluids to be immiscible and approach equilibrium quickly whilst ensuring 
numerical stability and positive shear viscosity. This amounts to keeping p as high as possible, 
r close to 1/2, and g as large as allowed by the onset of numerical instabilities which set in 
when the forcing term is too large. A large g allows for the early time transient, dominated 
by diffusion, to be of short duration. Finally, seeking the diffusive regime means looking at 
very early times, which is attained for large values of To. Conversely, the hydrodynamic inertial 
behaviour requires as small values of To as possible. 

In Table ^ we present the parameters selected in this study, along with the measured surface 
tension. We also include the length and time scales associated with them, which are used to 
compute dimensionless lengths and times in the model. 



Parameter set 


P 


r 


9 


<r(p,T,g) 


Lo(p,T,g) 


T (p,r,g) 


I 


0.8 


2.000 


0.06 


0.002059 


97.1 


18870 


II 


0.8 


1.500 


0.06 


0.004777 


18.6 


1038.8 


III 


0.8 


1.000 


0.06 


0.010292 


2.16 


28.0 


IV 


0.8 


0.625 


0.05 


0.017458 


0.0796 


0.152 



Table 1: Model parameters studied, including the surface tension a(p,r,g) measured for a planar inter- 
face on a 4 x 4 x 128 lattice, and the characteristic length Lq and time To for each parameter set. The 
existence of the latter two is based on the validity of the dynamical scaling hypothesis, and that diffusive 
currents are negligible with respect to hydrodynamic currents and capillary forces. 

The initial condition used for all the simulations was a thorough mixture of the two phases, 
with randomly distributed fluctuations. To realise this, each velocity direction k at each lattice 
site was populated with one density p% (x, t) = m a n^(x, t) for each species a = R, B as a white- 
noise, pseudo-random floating point number between 0.0 and 0.8, where m a are the particle 
masses, all set to unity. Note that the density p in Tabled is defined as the lattice average 

p=(p*(x,t) + p B (x,t)), (37) 

where p a (x,t) = J2k Pk( x '^)' an< ^ due ^° * ne critical composition we use, amounts to the maxi- 
mum value of either of the summands. 

Lattice sizes used were 128 3 and then 256 3 to check for finite size effects. Simulations for 
128 3 systems were run for 700 or 1400 time steps, and for 200 or 250 time steps for 256 3 systems, 
depending on the parameter set. Following Kendon et al.'s prescription to keep finite size effects 
at bay [Jjj, we neglected domain sizes larger than a quarter of the lattice side size. There is no 
reason a priori to choose this particular threshold. As we shall see, this allows the generation 
of a domain size range broad enough for data acquisition; furthermore, finite size effects were 
quantified by using the two aforementioned lattice sizes. 

Surface tension was measured on 4 x 4 x 128 and 16 x 16 x 128 lattices, allowing plenty 
of room along the non-isotropic direction z for the fluid's physical quantitites to achieve values 
characteristic of the bulk before being affected by the second interface with periodic boundary 
conditions imposed. We found that the surface tension did not vary by more than 1% when the 
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length along the z direction was doubled, which is the only direction where we would expect 
any variation as translational symmetry is broken. 

To compute the average domain size, Eq. (j2J), we perform discrete Fourier transforms. 
The sampling theorem |63| warns us to ensure that our fluid mixture does not exhibit spatial 
frequencies larger than the Nyquist critical frequency f c , defined as half the sampling frequency. 
This not being the case, the power spectrum in the interval [0, f c ] is altered by frequencies 
larger than f c as a result of aliasing. Because the sampling frequency on the lattice is one, the 
maximum frequency any relevant quantity of our fluid mixture is allowed to have according to 
the sampling theorem is 1/2, i.e., of wavelength two. This means that any spatial variation is 
bound to happen between two contiguous lattice sites, which is something we already knew: 
the resolution of the lattice is finite and dictated by the lattice size. We used the FFT routine 
rlft3() for real, 3D data sets [63] . 

Calculation of the reduced time t requires an assessment of Tint, ^int serves to redefine the 
time such that the domains have zero size at the time origin, which is not the case in the actual 
simulations. Depending on the regime reached by the parameter set employed, domains may 
start to grow immediately after time step zero, completely avoiding the diffusive stage. 

We assess T n t in the following way. We first compute the intersection with the abscissae of 
a linear fit interpolating all data starting after the initial purely diffusive transient is completed, 
that is, for which interfaces are thin enough and L(T) just starts to grow. The intercept is used 
as an initial guess for a\ in a Levenberg-Marquardt non-linear least-squares fit of the form 

y = a (x - ai ) a2 . (38) 

Once Ti n t is computed, and the data sets are normalised by Lq and To, hence becoming (t, I) 
data pairs, we perform fits to the function (|38[) to determine the growth exponent 02- Initial 
guesses for the fitting coefficients are a(j = 1.0, a? = 0.1, and a® = 1.0. The tolerance for these 
fits was set to 10 -5 , this being defined as the unsigned increment of x 2 between two consecutive 
iterations, divided by the number of degrees of freedom. 

Uncertainties in parameters are also taken care of. Because standard errors {A^S 1 } are 
incurred in the structure function spherical averaging these transmit down to L and I, 
and to t through the determination of T; nt . In this study, however, errors in the abcissae are 
disregarded as they do not depend on time, and therefore represent equal weights for data points 
in the least-squares functional to minimise. 

We performed the simulations using a number of processors ranging from 32 to 128 on a Cray 
T3E-1200E and on SGI Origin2000 and Origin3800 supercomputers. The code is an implemen- 
tation in Fortran90 using the Message Passing Interface (MPI) as parallelisation protocol, and 
it shows scaling with the number of processors between 50% to 90% of linearity on the Cray 
T3E platform up to 64 PEs, and better behaviour on SGI Origin platforms. CPU times 
used up to run a 128 3 lattice for 1400 time steps, or a 256 3 lattice for 250 time steps, ranged 
between 3 to 6 hours per parallel process. 

An important issue in dealing with the lattice sizes employed here is to have access to massive 
disk storage. For our largest lattices, 1.9 Gbytes of measurements were dumped onto disk at 
each measurement time step. A lattice of 256 3 sites run for 700 time steps, measuring every 25, 
requires 40 Gbytes to store the order parameter, the density fields for each phase, momenta, and 
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checkpoint files, the latter being needed if we wish to restart the simulation at the point where 
it stops. To that we need to add some additional working space for converting the dumped 
binary data into machine-portable XDR format |64| . For this work we required 200 Gbytes on 
disk, plus tape storage to free up space when required. XDR files were visualised using the 
commercial package AVS |65| . 

It is worth noting that our results did not undergo a process of lattice size reduction, in the 
sense of averaging over nearest-neighbouring sites in order to deal with limited computational 
resources, as was done in previous studies on 3D spinodal decomposition [201 EH Hence, we 
benefitted from measuring and visualising all data output from our simulations. Current limita- 
tions in computing resources prevented us from simulating lattices of 512 3 or 1024 3 sizes, which 
would otherwise be desirable in order to decrease the fluid's minimum Knudsen number, help- 
ful in reaching for the thermohydrodynamic limit as a multiscale Chapman-Enskog expansion 
procedure shows. However, this situation is bound to change soon with the advent of terascale 
computing capabilities. 

5.1 Growth exponents 

Figure^shows the average domain size in lattice units as obtained straight from the simulations, 
for all parameter sets (c/. Table Reynolds numbers achieved for each of these are Re = 0.18, 
0.49, 2.7, and 37. For parameter set I, we can see that after a transient during which there is a 
rapid mass convection to nearest neighbours, domain growth flattens out and starts growing at 
about T = 400. We will look at this in further detail; for now it can be seen that the breadth 
of the plateau decreases with the Reynolds number. Finally in Fig. [2] we show the same curves 
after rescaled to Lq and To, in reduced units. 

Fits to the model y = ao(x — ai) a2 for Fig. [2]are given in Tabled and they proved to be quite 
sensitive to the number of points fitted. Domain growth shows an increasing segregation speed, 
t c '• , t°- 593 5 t - 623 , and i 0,717 , with increasing Reynolds number. These data sets correspond 
to characteristic lengths and times in the ranges 0.0796 < Lq < 97.1 and 0.152 < Tq < 18870. 
These contain the values for which Kendon et al. [5] observed a viscous linear exponent, Lq = 5.9 
and Tq = 71. This, therefore, invalidates the universality of the dynamical scaling hypothesis. 

By looking at Fig. [2] from a grazing angle one can easily see that a simple, algebraic inter- 
polating curve is not obtainable here. Kendon et al. [HI EE1 and Pagonabarraga et al. [Ml EZJ 
used a method to improve this curve. They left T; nt as an adjustable fitting parameter such 
that there is a reasonable collapse onto a simple, single algebraic curve for all parameter sets 
simulated; from this they obtained a window of T| n t in which collapse is reasonable. Then they 
checked whether the different values for T[ n t from each individual parameter set lay within such 
a window. Quoting Kendon et al. (cf. Section 9.3 in Ref. [2S])) "although this [procedure] is 
capable of falsifying the scaling hypothesis [...], its non-falsification [...] may not represent per- 
suasive proof that the scaling is true." We adhere to this comment and prefer not to manipulate 
the data sets in such a way. 
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Figure 1: Evolution of the average domain size for parameter sets I, II, III, and IV (c/. Table ^) 
with the time step. Error bars are included and represent the uncertainty transmitted from the 
standard error of the structure function spherical average. Lattice size is 128 3 . All quantities 
are reported in lattice units. 
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Figure 2: Log-log plot of reduced length versus reduced time for the 128 3 -lattice data sets. 
Error bars are included. The four data sets correspond to parameter sets I, II, III, and IV (c/. 
Table from left to right. Viewed from a grazing angle, one can see that a simple, algebraic 
interpolating curve is not truly obtainable here. The first few points of each set correspond to 
diffusive, zero-growth stages. The units on both axes are dimensionless. 
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Parameter set 


«0 




0.2 


X /ndt 


I 


U.b44±(J.U14 


— z x 1U ±0.002 


n C/ic i r\ r\~\ A 

U.545±U.U14 


U.4b 


T T 

1 1 


n QO/i-Ln nn/i 

u . y L 4fc zt u . U U 4 


fi v in - 6_i_n nn7 


n f;n7-i-n nnfi 






0.922±0.003 


-2 x 1(T 5 ±0.007 


0.593±0.007 


0.48 


III 


1.248±0.031 


-0.007±0.100 


0.650±0.007 


2.71 




1.362±0.010 


-1 x 10" 4 ±0.03 


0.623±0.002 


0.68 


IV 


0.941±0.019 


0.01±3.9 


0.743±0.002 


0.10 




1.139±0.017 


-0.01±3.6 


0.717±0.002 


0.14 



Table 2: Levenberg-Marquardt nonlinear least-squares fits of I vs t data to model (|38(l . for each param- 
eter set attempted. The first line for each set belongs to 128 3 data, the second line to 256 3 data, the 
latter being unavailable for set I. Fitting parameters are given, plus the weighted sum of squared residuals 
(x 2 ) divided by the fit's number of degrees of freedom. Weights are the inverses of squared uncertainties. 
Note that x 2 / n df, also referred to as the variance of residuals, is expected to approach unity. Values 
larger than 1.0 may be due to an insufficient number of data points, data errors not normally distributed, 
or an incorrect model function. Values smaller than 1.0 may be the result of too large error bars, or too 
general a model function. 



5.2 Structure function 

For parameter set I ( cf. Table ^) we show in Fig. 01 a family of spherically-averaged structure 
functions versus wavenumbers, corresponding to time steps 200, 400, 600, 800, 1000, 1200, and 
1400, from right to left. Just as in scattering cross-section measurements |17| . we observe the 
peaks to grow and approach small wavenumbers as time evolves. In Figs. 21 and [5] we show the 
same family of curves using time steps as absissae and wavenumbers as parameters. Regions of 
linear growth with time on such a logarithmic scale indicate that a diffusive process is dominating 
the dynamics. In fact, an exponential time growth for the structure function shortly after 
the quench below the spinodal curve was predicted from the linearised Cahn-Hilliard Model-B 
equations without noise which although incorporating order-parameter conservation, does 
not include hydrodynamics. This Cahn-Hilliard equation might be applicable to regimes in our 
fluid where hydrodynamic effects were unimportant, as in the initial stages. Assuming linear 
perturbations <p' to the order parameter, Cahn predicted that for fluctuations of small amplitude 
and long wavelength there is an instability of the shape 

S(k, t) = S(k, 0) e - 2 ^* (39) 
for k < k c , where k c depends on the diffusion constant. Here, t is the time, to{k) < 0, and 

2 

S(k,t) oc ( </>' k (£) ), the brackets denoting averaging in reciprocal space over a shell of radius k. 

Exponential growth occurs in our simulations as can be seen from Figs. 0and[21for about the 
first 350 time steps for most of the wavenumbers, indicating its transient character. The plateau 
of Fig. ^ Set I, lasts during the first 400 time steps, and we can see, Fig. that up to 400 
time steps the peak in the structure factor varies in height and very little in wavenumber, and 
is located at 0.491 (lattice units). This leads us to think that at these early stages the dynamics 
is mainly making walls thinner while average domain sizes barely change. In addition, visual 
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inspection of the order parameter confirms the latter and suggests that hydrodynamic currents 
are weak, leaving diffusion as the mechanism leading the phase segregation process. When we 
check the structure function temporal evolution, Figs. 0] and for the curves at and around 
k = 0.491, we see that up to exactly 400 time steps do they show exponential growth, as the 
Cahn-Hilliard Model B predicts for a diffusive scenario. Also, exponential growth does not hold 
for all wavenumbers, but only for those smaller than about 0.7, in agreement with the existence 
of an upper cutoff for the validity of Eq. (|39|) . predicted from Model B. 

However, not all the wavenumbers follow Model B's predictions, namely, that exponential 
growth is a transient and occurs up to a threshold wavenumber. In fact, exponential growth holds 
for all the time steps of the simulation for the larger length scales (wavenumbers up to about 
0.245), suggesting that diffusion never ceases to dominate their dynamics. Also, exponential 
growth is seen for very small domain sizes (wavenumbers larger than 0.736) for time steps well 
advanced in the coarsening dynamics, after 600 time steps. These wavenumbers are close to 
and above the expected Model-B upper cutoff for exponential growth, set by the change in 
slope from positive to negative in Fig. [5J These departures from Model B's predictions hold 
nonetheless for domain sizes far from the first moment of the structure factor, which is close 
to its peak and is our average domain size measure. It would be desirable in future works to 
investigate diffusional processes at k < 0.245 for all of the simulation time, and k > 0.736 at late 
times: according to the Cahn-Hilliard linearised Model B, for these cases diffusion is negligible 
or forbidden, respectively. 

Analogous behaviour to Fig. 01 is exhibited for parameter sets II, III, and IV (c/. Table 
in Figs. EJCl and|Hl respectively. For the last two time slices taken in Fig. |5J the peaks seem no 
longer to drift to the left, as a result of finite size effects (arrest of domain growth). Regarding 
regions of exponential growth with time, the three data sets confirm Eq. (|39|) . with an upper 
bound for k. 

Figure |H1 shows the collapse (matching) of the structure functions corresponding to parameter 
set III (c/. Tabled, for a 128 3 lattice size and time steps from 450 to 700, when they are scaled 
by Eq. Q, the abscissae being rescaled by a factor of (27r) _1 , and the ordinates by the peak's 
maximum. Earlier times are represented in Fig. [HI by empty symbols, and later times by filled 
symbols. There is good collapse, and, therefore, scaling according to the scaling hypothesis, in 
the region from q = 0.4 to about q ~ 3, where q = kL is dimensionless. The middle of the region 
1 < q < 2 follows a q~ 9 behaviour, in accordance with Tomita's prediction of an exponent —6 
or more negative |69j . 

Close to q ~ 3 we observe the presence of a shoulder, as has been reported in experiments 
|67| and numerical simulations (541 168j . Most strikingly, the shape of our large-g tail is very 
reminiscent of that of Fig. 4 in j^I] and that of Fig. 3 in [HE]: (1) there is still a time dependence 
indicating that interfaces have not yet been fully resolved (we are probing the smallest scales 
where diffusion still exists and £/L is not small enough); and (2) the tail decreases with an 
exponent which is in fact more negative than that of the Porod tail, Eq. ©, despite what these 
authors jSUISH] claim. 

For q < 0.4, data points do not seem to collapse onto the same curve of those for q > 0.4. 
This is similar to, but with more data than, the results of Koga and Kawasaki [HSJ. Our results 
show an exponent growing with time: the slope of a line (not shown) joining the first two 
empty circles (T = 450) is 1.61, while the slope of a line (not shown) joining the last two filled 
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downward triangles (T = 700) is 2.12. This resembles the temporal growth cited by Appert et 
al. |54j on the results of Alexander et al. |35j : nonetheless, we consider the amount of data 
in the latter insufficient to draw firm conclusions. Given that the points at T = 700 are closer 
to the asymptotic regime, we take such a slope as our best approximation to the asymptotic 
regime. 

In the small-g region, Yeung [H^ predicted a q behaviour for the asymptotic limit (L —* oo, 
or at late times). Additionally, at earlier stages, a term proportional to L~ 2 q 2 caused by thermal 
noise would also come into play. Now, Appert et aUs estimate [3] applies well for our results: 
such a quadratic term is less dominant than the quartic one only for q > 0.4, given that the 
largest value of L(T) for which there are no finite size effects is also about 25. This happens to 
be the region where we find the q 2 <-> q 4 crossover. 

Fig. El] shows similar curves for parameter set IV (c/. Table where only time steps 450 
to 675 are displayed and we have also normalised the curve such that the peak is located at 
(1,1). We have again neglected early time steps because of poor collapse. A fit to the tail in 
2 < q < 10 gives g _3 - 65 ; close to being a Porod's law. It is when we probe the finest length 
scales, at q « 10 that it ceases to apply, due to lattice discretisation effects. 

The behaviour at intermediate wavenumbers is between q~ 8 and q~ 7 , again in agreement 
with Tomita's theory jSH], and close to q~ 7 as computed using a dissipative particle dynamics 
method by Jury et al. an d a lattice-gas automaton by Love et al. [Zj. 

For small momenta (large domains) we found a behaviour close to q 3 , in agreement with the 
numerical results of Love et al. [J] and in disagreement with Yeung's predictions |61| . 

The most notable difference between Figs. 151 and 1TU1 is the behaviour above q ~ 1.5. Figure 
1101 shows a neat Porod tail, which bends down dramatically for q > 10, whereas Fig. |8j shows 
either a poor Porod tail in the region 3 < q < 5, or a minute one in the region 1.5 < q < 3. A 
condition assumed in the derivation of Porod's law JH] is that the sampling length r satisfies 
£ -C r -C L, which in wavenumbers means 



By 'eyeball' inspection of the system's order parameter we found that interface widths naturally 
shrink with an increasing number of time steps, going from about 5 or 6 lattice unit spacings 
at 200 time steps down to about 3 at 675 time steps, regardless of the data set, III or IV (c/. 
Table Simulations for a 256 3 lattice size revealed similar widths, and snapshots of the order 
parameter at 200 and 700 time steps are shown in Figs. llll and lT2l With these widths in mind, 
assuming domain sizes of a quarter of the lattice side length (the threshold imposed by our 
prescription for eliminating finite size effects), and a 128 3 lattice, the condition (|40[) becomes 



which contains our large-g region. Despite this, we do not observe a Porod tail for data set III, 
or, as in data set IV, the tail obtained is only close to being a Porod tail. This is in agreement 
with the fact Eq. |lU]is necessary but not sufficient for a Porod tail to hold. 

Finally, it is worth noting in Fig. El the existence of nested domains and droplets much 
smaller than the average domain size. 



l/L < k < l/£. 



(40) 



1«?«10, 



(41) 
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Figure 3: Spherically averaged structure function versus wavenumber, for parameter set I (c/. 
Table^). 128 3 lattice. Error bars represent the standard error of the structure function spherical 
average. Time slices shown are time step 200, 400, 600, 800, 1000, 1200, and 1400 from right to 
left. All quantities are reported in lattice units. 
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Figure 4: Evolution of the spherically averaged structure function with the time step for param- 
eter set I and a 128 3 lattice, on a logarithmic scale. When observed along the ordinate T = 200, 
the curves correspond to wavenumbers k = 0.147, 0.196, 0.245, 0.295, 0.344, 0.393, and 0.442 
from bottom to top, respectively. Error bars represent the standard error of the structure func- 
tion spherical average. Regions of linear growth are those for which the exponential behaviour 
Eq. (|39|) holds. For wavenumbers up to 0.2 exponential and therefore diffusive behaviour is seen 
for all the simulation time. For larger wavenumbers (and hence smaller domain sizes) diffusion 
occurs as a transient. All quantities are reported in lattice units. 
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Figure 5: Similar to FigEJ but the curves correspond to wavenumbers k = 0.491, 0.540, 0.589, 
0.638, 0.687, 0.736, 0.785, and 0.834 from top to bottom, respectively. We can see that linear 
growth ceases to hold for wavenumbers larger than about 0.736, in accordance with existence of 
an upper cutoff for the validity of Eq. (j3*9"j) . All quantities are reported in lattice units. 
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Figure 6: Spherically averaged structure function versus wavenumber, for parameter set II (c/. 
Table^). 128 3 lattice. Error bars represent the standard error of the structure function spherical 
average. Time slices shown are time step 200, 400, 600, 800, 1000, 1200, and 1400 from right to 
left. All quantities are reported in lattice units. 
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Figure 7: Spherically averaged structure function for parameter set III (cf. Table^). 128 3 lattice. 
Error bars represent the standard error of the structure function spherical average. Time slices 
shown are time step 100, 200, 300, 400, 500, 600, and 700 from right to left. All quantities are 
reported in lattice units. 
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Figure 8: Scaled spherically averaged structure function for parameter set III (c/. Table as 
defined by Eq. (jHJ). Lattice size is 128 3 . Time steps are as shown in the legend. Earlier times 
correspond to the empty symbols; later times to the filled symbols. Error bars are smaller than 
the size of the symbols. Straight lines serve as slope guides to the reader only, and represent 
power laws q 2 , q 4 , q~ 9 , and Porod's law q~ 4 , from left to right, respectively, with q = kL. All 
quantities are reported in lattice units. 
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Figure 9: Spherically averaged structure function versus wavenumber, for parameter set IV (c/. 
Table^). 128 3 lattice. Error bars represent the standard error of the structure function spherical 
average. Time steps shown are 100, 200, 300, 400, 500, 600, and 675, from right to left. All 
quantities are reported in lattice units. 
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Figure 10: Scaled spherically averaged structure function for parameter set IV (cf. Table 
as denned by Eq. (J8J). Lattice size is 128 3 . Time steps shown are from 450 up to 675, every 
25. Earlier times correspond to the innermost lines; later times to the outermost lines. Error 
bars are smaller than the size of the symbols, except for the two leftmost, detached data sets, 
for which they are slightly larger. Straight lines serve as slope guides to the reader only, and 



represent power laws q , q , q , q , and the fit to the large-g tail, q 
respectively, with q = kL. All quantities are reported in lattice units. 
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Figure 11: Order parameter (p — p ) snapshot at time step 200 for parameter set IV (c/. Table 
P). We show a slab of the 256 3 lattice simulated. 
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6 Numerical stability of our lattice-Boltzmann algorithm 



As is well known, owing to the lack of an H-theorem, an approach to equilibrium is not guar- 
anteed in all lattice-Boltzmann models to date; recent theoretical developments to address and 
solve this have been made |421 1431 IH] . For single-phase lattice-Boltzmann models, equilibrium 
states are well defined in the collision term; if the automaton does relax to these, the perti- 
nent macroscopic momentum (and sometimes energy) balance equations are reproduced in the 
low-Knudsen number limit. Interacting, multicomponent lattice-Boltzmann models exhibit the 
same situation in the bulk of pure fluid regions where inter component interactions are negligible. 
For regions where they are not, there is not even a well-established thermohydrodynamic theory 
which could provide equilibria to which the automaton could relax to, or with which to com- 
pare the stationary state to which it can evolve. Whether dealing with a single or multiphase 
system, lattice-BGK stationary regimes ought to be treated with caution and contrasted with 
experiment. 

Numerical instabilities are the reflection of the lack of an H-theorem, which is a direct 
consequence of space and time discretisation on the BGK-Boltzmann equation and the freedom 
in the choice of the equilibrium distribution function [H] . These instabilities can be defined as 
follows. As is generally the case for a finite difference method with a single relaxation parameter, 
such as our lattice-BGK model for a zero phase-coupling constant jlj, linear stability occurs 
within a finite interval of such a parameter. If multicomponent interactions are introduced, 
additional parameters may influence the stability: density, intercomponent coupling strength, 
and even composition. The mechanism is simple: certain choices of parameters can turn the 
lattice-BGK collision term positive (therefore increasing the mass density) for long enough to 
generate floating-point numbers larger than the largest the machine can deal with, hence causing 
an overflow signal. Numerical instabilities are defined in this work as the generation of such 
floating-point numbers. We consider it crucial to be able to map out regions in the model's 
parameter space leading to unstable configurations, and to report them alongside any lattice- 
BGK simulations. 

Using the same initial conditions as explained in Section |S1 we found our algorithm to be 
unstable for regimes with the smallest length and time scales, Lq and To, which coincide with 
those of the largest Mach numbers. In Table 131 we show some of the paramaters leading to 
numerical instability. The dependence of the surface tension on the model parameters, as given 
in Section |IJ should be taken into consideration as a guide to steering through the parameter 
space. Note that all values of At included are larger than that for parameter set IV. 

We then investigated the nature of our instabilities, as others have done. The group of Cates 
found troublesome numerical instabilities with their free-energy based, lattice-BGK model in 
3D in regions in which quiescent binary portions of fluid go into a checkerboard state |27j . 
They reported that their model is unconditionally unstable |26| . Nonetheless, by improving the 
way gradients were treated numerically they were able to considerably reduce this unphysical 
behaviour. For our model, we looked at the time evolution of the quantity 

9{t) = max{|Qf(x,i)| , Vx,VA;,Va} , (42) 
for parameters {p = 0.3, g = 0.06, r = 0.5125}, where the collision term, f2C™, is defined in Eq. 
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Table 3: Model parameters leading to numerical instability, including the surface tension a(p, t, g) 
generated some time steps before the instability sets in, and the associated characteristic time. The 
lattice used was 4 x 4 x 128, and the instability sets in before 4000ts. 

(|25|). We also monitored the maximum and average values of the fluid mixture's speed, t£ max 
and u, respectively, on the lattice. We show these quantities for a 32 3 lattice in Figs. IT3l ITU and 
H31 We see how 6 reverses its decreasing trend in a few time steps; after that, it blows up at 
T = 52 time steps. We only show data up to T = 49, as 9(T = 51) « 10 111 . u max blows up in 
similar style: at T = 50, "Umax = 

7498 and d « 10 21 ; at T = 51, K max has exceeded the maximum 
floating-point value that the computer can deal with, and overflow signals are generated. This 
indicates that at the time steps immediately prior to the onset of the instability the lattice gets 
more and more populated with increasing speeds until in two or three time steps they grow by 
ten or more orders of magnitude. That the population of lattice sites with rapidly increasing 
speeds over time is small compared to the lattice volume can be concluded from contrasting the 
time variation in the standard error ( "one sigma" ) of u to the time variation of u, Fig. 1151 The 
same parameter set run on a 128 3 lattice seems to make the instability set in much quicker, as it 
occurs during the first 10 time steps. As a final check, we ran a 128 3 lattice with parameter set 
I (cf. Table ^) for 20000 time steps and found no instabilities. The time evolution of 6, u and 
u is shown in Figs. Cxi and ^1 respectively. We conclude that the occurrence of instabilities 
only depends on the set of parameters used, regardless of the number of time steps simulated. 

7 Conclusions 

We have presented a quantitative study of the phase separation dynamics in three dimensions for 
critical (50:50) fluid mixtures (spinodal decomposition) for a modified Shan-Chen lattice-BGK 
model of multicomponent, isothermal immiscible fluids. 

We found that after a brief diffusional transient in which interconnected regions of fluid 
species embedded into one another are formed, the average size of such regions grows with time 
as I oc i 7 , where 7 ~ 2/3. The trend is for the value to increase in the range 0.545 ± 0.014 < 
7 < 0.717 ±0.002 as the Reynolds number increases. This increase is consistent with a crossover 
from I oc i 1 / 3 diffusive behaviour to hydrodynamic viscous growth I oc t predicted by the Cahn- 
Hilliard Model H. Owing to the significant amount of diffusion present at low Reynolds number, 
we do not consider our results to be indicative of a genuinely hydrodynamic inertial t 2 ^ 3 regime. 

We observed exponential growth in the time dependence of the structure function for wavenum- 
bers up to a threshold value, in qualitative agreement with predictions from the linearised Cahn- 
Hilliard Model B. For small wavenumbers, such an exponential growth is seen at all simulation 
times, whereas it is only an initial transient for larger wavenumbers. These departures from 
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Figure 13: Evolution of the collision term maximum absolute value, 9, Eq. 1)42)1 . with the time 
step on a 32 3 lattice for parameters {p = 0.3, g = 0.06, r = 0.5125}. The interpolating curve 
serves as a guide to the eye only. All quantities are in lattice units. 
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Figure 14: Evolution of the maximum speed, « max , with the time step on a 32 3 lattice for 
parameters {p = 0.3, g = 0.06, r = 0.5125}. The interpolating curve serves as a guide to the 
eye only. All quantities are in lattice units. 
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Figure 15: Evolution of the speed average, u, with the time step on a 32 3 lattice for parameters 
{p = 0.3, g = 0.06, r = 0.5125}. Error bars represent the standard error of the average ("one 
sigma"). All quantities are in lattice units. 
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Figure 16: Evolution of the collision term maximum absolute value, 9, Eg. 1)42(1 . with the time 
step on a 128 3 lattice for parameter set I (cf. Table We can see a decreasing trend for most 
of the simulation, which accentuates after time step 10000. The interpolating curve serves as a 
guide to the eye only. All quantities are in lattice units. 
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Figure 17: Behaviour of the maximum speed, u max , with the time step for a 128 3 lattice with 
parameter set I (cf. Table ^). It shows an overall decreasing trend. The interpolating curve 
serves as a guide to the eye only. All quantities are in lattice units. 
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Figure 18: Evolution of the speed average, u, with the time step on a 128 3 lattice for parameter 
set I (cf. Table ^). We can see a decreasing trend of the average and its error. All quantities 
are in lattice units. 
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Model B predictions are for wavenumbers far from the one characterising the average domain 
size. A natural continuation of this work would be to investigate the nature of diffusion currents 
for these cases. 

We have found very good agreement with the dynamical scaling hypothesis in the form 
of a neat collapse of the structure function curves for Re = 2.7 and Re = 37 when they are 
appropriately scaled according to Eq. l|K|). This collapse holds roughly for the second half of the 
simulation time, as diffusional transients act during the first. By looking at order parameter 
snapshots we observed the formation of nested domains and smaller droplets for the largest 
Reynolds numbers achieved, as Wagner and Yeomans also found |1U| . However, unlike them, in 
our case these are transients rather than a result of length scales growing at different speeds, as 
poor collapse of the scaling functions would then occur due to breakdown of scale invariance. 

Yeung predicted a q 2 behaviour at the small-g end of the spectrum as the result of thermal 
effects at pre-asymptotic stages [H^. Love et al. [?] conjectured that a q 2 behaviour, and a 
crossover to q 4 , could be caused by (a) lattice-gas noise, or (b) a poor scaling collapse, and 
that their i 2 / 3 domain growth, instead of t, might be justified by the former. Appert et al. 
|54j ascribed the q 2 behaviour and the crossover to not having reached the asymptotic limit, 
L — > oo (poor scaling collapse again). Our noiseless model reproduced the q 2 «-> q 4 crossover at 
Re = 2.7 and did not at Re = 37, for which there is better scaling collapse, and also produced 
a 2/3 domain growth (crossover) exponent. All this leads us to conclude that noise may not 
play as important a role as the lack of scaling collapse in explaining the q 2 <-> q 4 crossover, 
and is definitely not a requirement for the reproduction of a 2/3 domain growth exponent. A 
q 2 behaviour is the only one experimentally observed by Kubota et al. .67 for a mixture of 
isobutyric acid and water; they cite surface tension effects, measurement difficulties, multiple 
light scattering, and even specificity to the mixture's molecular weight as reasons for not seeing 
a q 4 behaviour, and definitely discard thermal noise. Not surprisingly, in his prediction Yeung 
assumed a diffusive domain growth exponent of 1/3, which is rather seen in quenches of polymer 
mixtures and metal alloys. 

In the case Re = 37, the spectrum shows a q s behaviour in the small-g limit, in disagree- 
ment with Yeung's prediction. In fact, his analysis is based on a Cahn-Hilliard model without 
hydrodynamics. 

The numerical instabilities seen in our runs are caused by large speeds turning the equi- 
librium distribution negative for long enough to incur floating-point overflows. This happens 
for characteristic times (c/. Table ^) below To = 0.0172, and the population of lattice sites 
undergoing such a burst in the fluid's macroscopic speed is small compared to the lattice vol- 
ume. We found no evidence that an initially stable regime becomes unstable at later times, as 
typically happens in relaxational models (such as is our model for g a a = 0). This is in stark 
contrast with the findings of Kendon et al. [HH and Cates et al. j^Tj in their spinodal decom- 
position studies using a free-energy based, lattice-BGK model, who reported their algorithm to 
be unconditionally unstable. 

A search for a crossover to growth laws other than i 2 / 3 at Reynolds numbers higher than 
Re = 37 faces two major problems: (a) the triggering of numerical instabilities due to large 
inter-species coupling and smallness of relaxation time values; and (b) the approach to the 
compressible limit, whose macrodynamic behaviour for pure phases cannot, by construction, be 
correctly described by our method. On the other hand, there is still scope to achieve Reynolds 
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numbers smaller than Re = 0.18 in search of the end of the crossover to i 1 / 3 . Closeness to the 
miscibility threshold may make this endeavour difficult, as it is reached for characteristic times 
ca. T = 1.43 x 10 8 . 

Our results clearly challenge the claim that a domain growth linear with time is universal 
for all models of phase separating fluids sharing similar values of Lq and To since we obtained 
excellent collapse of scaled structure functions yet our domain growth exponents are in the 
crossover region between diffusive and hydro dynamic viscous regimes. 

The properties of this binary immiscible fluid model are important for underpinning the 
more complex domain growth observable in ternary amphiphilic (oil/water/surfactant) fluids 
we expect to report in forthcoming publications. 
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